###NMDS
setwd("D:/入侵/数据/分析2021/3 NMDS")
###NMDS
species<-read.csv("species.csv",header=T)  
rownames(species)<-species[,1]  #Set the first column to the row name
species<-species[,-1] #delete the first column
species_no_com<-species[-1,] #delete the first row（Community）

species<-species_no_com[,-1] #delete the original second column（trophic）
dat=as.matrix(species)
dat=matrix(as.numeric(dat),nrow=nrow(dat),dimnames = dimnames(species))

dat_t<-t(dat)
map<-read.csv("mapping.csv",header=T)
sol <- metaMDS(dat_t)

sim=simper(dat_t,map$Community)

sim_result_OP_RP=data.frame(sim$OP_RP$ord,sim$OP_RP$cusum)
write.csv(sim_result_OP_RP,file = "sim_result_OP_RP.csv")

sim_result_RP_IS=data.frame(sim$RP_IS$ord,sim$RP_IS$cusum)
write.csv(sim_result_RP_IS,file = "sim_result_RP_IS.csv")

sim_result_RP_IS=data.frame(sim$OP_IS$ord,sim$OP_IS$cusum)
write.csv(sim_result_RP_IS,file = "sim_result_OP_IS.csv")

sink("simper result.txt") 
print("OP_RP")
print(sim$OP_RP$overall)
print("OP_IP")
print(sim$OP_IP$overall)
print("OP_IS")
print(sim$OP_IS$overall)
print("OP_IPS")
print(sim$OP_IPS$overall)
print("RP_IP")
print(sim$RP_IP$overall)
print("RP_IS")
print(sim$RP_IS$overall)
print("RP_IPS")
print(sim$RP_IPS$overall)
print("IP_IS")
print(sim$IP_IS$overall)
print("IP_IPS")
print(sim$IP_IPS$overall)
print("IS_IPS")
print(sim$IS_IPS$overall)
sink()

NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],Site=map$Site, Community=map$Community)
NMDS$Community = factor(NMDS$Community, levels=c('OP','RP','IP','IS','IPS'))
p_NMDS=ggplot(data = NMDS, aes(MDS1,MDS2,group=Community, color = Community))+
  theme_customize+
  scale_color_manual(values=c("#44BBBB", "#f28a6a", "#8ea1ca", "#8CBB44", "#E6BD1A"))+
  theme(panel.grid=element_blank())+
  geom_point(aes(color = Community),size=3)+
  geom_text(label=NMDS$Site,size=3,hjust=1, vjust=-1)+
  theme(legend.key=element_rect(colour="white",fill="white",size=3,linetype="dashed"))+
  theme(legend.title=element_blank())+
  stat_ellipse(level = 0.95, show.legend = F) +
  annotate("text",x=min(NMDS$MDS1),y=min(NMDS$MDS2),hjust=0,vjust=0,label=paste("Stress:",sol$stress)) #6*5

anosim_all=anosim(dat_t,map$Community,permutations=999,distance="bray")


##ANOSIM
group_name <- unique(map$Community)
anosim_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {
  for (j in (i + 1):length(group_name)) {
    group_ij <- subset(map, Community %in% c(as.character(group_name[i]), as.character(group_name[j])))
    otu_ij <- dat_t[group_ij$Site, ]
    anosim_result_otu_ij <- anosim(otu_ij, group_ij$Community, permutations = 999, distance = 'bray')     
    
    if (anosim_result_otu_ij$signif <= 0.001) Sig <- '***'
    else if (anosim_result_otu_ij$signif <= 0.01) Sig <- '**'
    else if (anosim_result_otu_ij$signif <= 0.05) Sig <- '*'
    else Sig <- NA
    anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', 
                                                    anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif, Sig))
  }
}
#Outputs a table with R and P values
anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)
names(anosim_result_two) <- c('group', 'distance', 'R', 'P_value', 'Sig')
write.csv(anosim_result_two,file = "anosim_result_two.csv")

###trophic
trophic<-read.csv("trophic.csv",header=T)  
rownames(trophic)<-trophic[,1]  
trophic<-trophic[,-1] 
trophic_no_com<-trophic[-1,] 

trophic<-trophic_no_com[,-1] 
dat=as.matrix(trophic)
dat=matrix(as.numeric(dat),nrow=nrow(dat),dimnames = dimnames(trophic))

dat_t<-t(dat)
map<-read.csv("mapping.csv",header=T)
sol <- metaMDS(dat_t)

sim=simper(dat_t,map$Community)

sink("simper result for tro.txt")
print("OP_RP")
print(sim$OP_RP$overall)
print("OP_IP")
print(sim$OP_IP$overall)
print("OP_IS")
print(sim$OP_IS$overall)
print("OP_IPS")
print(sim$OP_IPS$overall)
print("RP_IP")
print(sim$RP_IP$overall)
print("RP_IS")
print(sim$RP_IS$overall)
print("RP_IPS")
print(sim$RP_IPS$overall)
print("IP_IS")
print(sim$IP_IS$overall)
print("IP_IPS")
print(sim$IP_IPS$overall)
print("IS_IPS")
print(sim$IS_IPS$overall)
sink()

NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],Site=map$Site, Community=map$Community)
NMDS$Community = factor(NMDS$Community, levels=c('OP','RP','IP','IS','IPS'))
library(ggplot2)
p_NMDS_tro=ggplot(data = NMDS, aes(MDS1,MDS2,group=Community, color = Community))+
  theme_customize+
  scale_color_manual(values=c("#44BBBB", "#f28a6a", "#8ea1ca", "#8CBB44", "#E6BD1A"))+
  theme(panel.grid=element_blank())+
  geom_point(aes(color = Community),size=3)+
  geom_text(label=NMDS$Site,size=3,hjust=1, vjust=-1)+
  theme(legend.key=element_rect(colour="white",fill="white",size=3,linetype="dashed"))+
  theme(legend.title=element_blank())+
  stat_ellipse(level = 0.95, show.legend = F) +
  annotate("text",x=min(NMDS$MDS1),y=min(NMDS$MDS2),hjust=0,vjust=0,label=paste("Stress:",sol$stress)) #6*5

anosim_all=anosim(dat_t,map$Community,permutations=999,distance="bray")


##ANOSIM
group_name <- unique(map$Community)

anosim_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {
  for (j in (i + 1):length(group_name)) {
    group_ij <- subset(map, Community %in% c(as.character(group_name[i]), as.character(group_name[j])))
    otu_ij <- dat_t[group_ij$Site, ]
    anosim_result_otu_ij <- anosim(otu_ij, group_ij$Community, permutations = 999, distance = 'bray')    
       
    if (anosim_result_otu_ij$signif <= 0.001) Sig <- '***'
    else if (anosim_result_otu_ij$signif <= 0.01) Sig <- '**'
    else if (anosim_result_otu_ij$signif <= 0.05) Sig <- '*'
    else Sig <- NA
    anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', 
                                                    anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif, Sig))
  }
}

anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)
names(anosim_result_two) <- c('group', 'distance', 'R', 'P_value', 'Sig')
write.csv(anosim_result_two,file = "anosim_result_two_tro.csv")